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1.  INTRODUCTION 


(1  2 

In  the  past  four  years  several  authors'  *  *  '  have  reported  Bulerlaa 
hydrodynamic  codes  for  the  computation  of  compressible  fluid  flows  in  which 
the  independent  variables  are  two  space  dimensions  and  time. 

iz's 

The  present  paper  describes  the  generalization  of  such  a  code'  7  to 
include  also  the  forces  due  to  strength  which  arise  within  the  material  to 
resist  shear  deformation.  In  choosing  the  precise  model  for  such  a  generali¬ 
zation  ,  it  is  useful  to  recall  that  Eulerlan  codes  have  proven  to  be  especially 
suited  to  the  description  of  flows  in  which  material  elements  undergo  extreme 
distortion.  It  is  in  precisely  this  area  of  treating  extreme  distortions 

that  Lagrangian  codes,  despite  their  success  in  treating  other  aspects  of 

(4  5) 

strength-dependent  deformation,'  encounter  serious  difficulties.  In 
order,  then,  to  fill  the  apparent  need  for  a  strength  code  which  is  capable 
of  the  simple  and  efficient  computation  of  flows  involving  large  deforma¬ 
tions,  the  present  code  development  was  initiated  in  1965 .  An  informal 
status  report  was  given, ^  and  the  code  has  subsequently  been  the  subject 
of  continuing  development  during  the  course  of  a  number  of  applications. 

A  basic  approximation  underlying  the  present  method  can  be  mentioned 
here.  Since  the  method  is  expected  to  have  its  primary  utility  in  the 
treatment  of  material  undergoing  large  deformations,  it  is  appropriate  to 
neglect  the  small  elastic  shear  strains  that  may  precede  the  onset  of  plastic 
deformation.  More  precisely,  for  a  given  increment  in  the  total  strain,  the 
de vie tor  part  of  the  strain  is  assumed  to  be  purely  plastic;  the  remaining 
dilatational  strain  is  fully  accounted  for  in  the  hydrostatic  equation  of 
state.  This  neglect  of  elastic  shear  strains  corresponds  to  a  so-called 
rigid-plastic  or  Levy- von  Mises  model  of  the  continuum,  to  be  discussed  as 
Section  V.  The  present  approximation  does  not  preclude  accounting  for  the 
dependence  of  the  material  strength  upon  the  thermodynamic  state  of  the 
material  or  upon  the  work  done  against  the  strength  forces  (work  hardening). 
These  effects  are  also  discussed  in  Section  5. 
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A  motivation  in  the  present  effort  was  the  desire  to  compute  the 
large  deformations  and  the  eventual  crater  size  occurring  in  the  hyper¬ 
velocity  impact,  and  the  results  02’  such  a  calculation  are  given  as  Ref.  It. 

For  completeness,  both  the  hydrodynamic  and  strength  aspects  of  the 
code  are  given  in  the  present  report,  although  the  hydrodynamic  part  is 
largely  a  description  of  the  code  which  was  reported  earlier. 
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2.  BASIC  EQUATIONS 


In  an  Eulerian  code  space  is  divided  into  fixed  cells  through  which  . 
the  fluid  moves.  To  arrive  at  expressions  for  the  rate  of  change  of  total 
mass,  momentum  and  energy  within  such  a  cell,  we  start  with  the  equations 
of  motion  in  the  form 

E  -  <’“!>  « 


»  Dt3  -  4  (•«>  <2> 

p  Dt  “  alj  (°ii  uj'  ^ 

Here  o^  is  the  stress  tensor,  which  can  he  regarded  as  the  sum  of  the  hydro¬ 
static  stress  -  &ijF  and  a  stress  deviator  tensor  s^,  i.e., 

au  *  Bu '  V  w 

c 

and  E  is  the  total  energy-  per  gram,  kinetic  plus  internal.  Tensor 
notation  is  implied,  so  that  repeated  indices  denote  summations. 

In  subsequent  sections  relations  will  be  given  for  determining  P 
from  the  equation  of  state  of  the  material  and  s^  from  the  material 
constitutive  relation. 

Expanding  the  convective  derivatives  in  nqs.  2-3,  Bf/Dt  =  df/dt  + 
u^  Bf/dx^,  then  adding  Eq.  1  times  u^  to  Eq.  2,  and  Eq.  .1  times  E  to  Eq.  3, 
and  collecting  terms,  gives 


it  (pV  "  ’ij  "  ^9ui  V 
it  (pE)  ‘  V  -  4  (puiE) 

For  the  developments  to  follow  it  is  desirable  to  replace  these 
differential  equations  by  related  integral  equations,  obtained  by 
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integrating  over  the  cell  volume  r,  and  then  converting  the  volume  integral 


*v 

of  divergences  to  surface  integrals  over  the  cell  surfaces. 
2%  and  3*  then  become 

Equations  1, 

&J“T  *iT  -  -I.  "Vt4" 

(5) 

h  U  "“A  ”  4  “Vi4*  -  4  Pttiu3ni4* 

(6) 

jt  Sr  ■  Is  VA4*  '  4  ’W' 
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3.  COMPUTATIONAL  METHOD 

3.1.  3MRpDUCTIQN_  -..DIVISION  .IOTOPHASES 

The  purpose  of  the  present  section  is  to  describe  the  computational 
method  by  which  the  flow  configuration  is  advanced  in  time.  This  is  done 
in  general  terms,  leaving  for  PART  THREE  the  specification  of  the  computer 
program  which  is  used  to  carry  out  the  actual  calculation. 

It  is  convenient  to  express  the  integral  conservation  relations, 

Eos.  5-7 ,  as  finite  difference  equations  over  the  time  step  At-  and  also  to 
decompose  the  total  stress  cr . <  into  its  deviator  and  hydrostatic  components, 
according  to  Eq.  4.  This  gives,  for  the  increments  of  total  mass  (m), 
momenta  (muj)  and  energy  (mE)  within  the  cell 

i 

Am  =  *  -  At  Jg  pu^ds  (8) 

A(m-^)  ■  -At  Js  Pn^ds  +  At  Sj^ds  -  At  Jg  (pu^Jn^s  (9) 

A(mE)  =  -At  Js  Pu^ds  +  At  J*g  s^u  n^s  -  At  Jg  (pu^J^ds  (10) 

Here  the  terms  on  the  right  are  divided  into  increments  due  to  the  pressure 
forces  on  the  cell  surface  (first  column),  those  due  to  the  stress  deviator 
forces  on  the  cell  surface  (second  column)  and  the  increments  (third  column) 
due  to  the  flux  of  mass,  momentum  or  energy  through  the  surface  of  the  cell. 
In  the  computation  these  three  types  of  contributions  are  accounted  for  in 
distinct  phases.  Specifically,  during  each  time  step  all  cells  are  first 
updated  for  pressure  effects  in  Phase  1,  the  effects  of  the  stress  deviators 
are  next  accounted  for  in  Phase  3  and  finally  the  transport  effects  are 
added  in  Phase  2.  In  the  discussion  that  follows  the  calculation  of  the 
terms  on  the  right  of  Eqs.  8-10  are  described  sequentially  starting  with 
Phase  1. 

Some  preliminary  definitions  will  be  helpful.  Superscript  N  on  a 
variable  refers  to  the  value  of  the  variable  at  the  beginning  of  the  time 
step  and  superscript  (N  +  l)  denotes  the  value  at  the  end.  In  this  first 
discussion  we  consider  a  typical  cell  in  the  interior  of  the  grid,  saving 
until  PART.'  3-6  a  discussion  of  the  special  conditions  which  are  required  at 
grid  boundaries  or  at  the  axis  of  symmetry.  For  a  typical  ce?l,  denoted 
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by  a  value  of  the  Index  K,  the  dependent  variables  for  that  cell  are 
written  P(K),  u(K),  v(k),  l(K),  M(k),  representing  respectively  the  pres¬ 
sure,  radial  and  axial  components  of  velocity,  the  specific  internal  energy 
and  the  mass  for  cell  K.  The  adjacent  cells  above,  below,  to  the  right  and 
left  of  K  will  be  designated  respectively  as  KA,  KB,  KR,  and  KL.  Here  the 
terms  above,  below,  right  and  left  refer  to  a  cross  section  view  of  the 
cells,  in  which  the  left  border  is  the  axis  of  symmetry  with  z  increasing 
upward,  see  Fig.  1.  Each  cell  is,  then,  the  torus  obtained  by  rotating  the 
rectangle  (since  hz  J  Ar  in  general)  about  the  axis  of  symmetry.  The  cells 
designated  by  KAL,  KAR,  LBL  and  KBn  will  be  referred  to  in  the  strain  rate 
calculation  discussed  in  Section  5. 


EDI 
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Fig.  1— Grid  layout  and  a  typical  cell 


3.2  PHASE  1  -  THE  EFFECTS  OF  PRESSURE 


3.2.1  Continuity  Equation,  Eq.  8 
No  contribution  in  Phase  1. 

3.2.2  Equations  of  Motion,  Eq.  9 

a.  Axial  Motion.  In  this  case,  referring  to  the  pressure  integral 
in  Eq.  9,  u^  *  v  and  n^  *  -1,  0,  +1  are  the  axial  components  of  the  unit 
normal  to  the  bottom,  sides  and  top  of  the  cell  respectively.  Equation 

9  therefore  gives  for  the  Phase  1  contribution 

A-^mv)  «*  Pb  «(rg  -  r^)  At  -  n(r g  -  r^)  At? 

where  rg  and  r^  are  the  radii  of  the  outer  and  inner  cell  surfaces, 
respectively,  and  P^,  P&  are  the  pressure  on  the  bottom  and  top  cell  faces 
These  pressures  are  obtained  from  the  initial  (time  N)  values  of  cell  pres 
sures  by  a  simple  average  of’  the  pressures  in  the  adjacent  cells-: 

PN(K)  +  PN(KB) 

'  b  =  2 

\ 

T  =  PN(K)  +  PN(KA)  . 

u 

Making  these  substitutions  gives,  for  the  Phase  1  increment  of  axial 
momentum 

Mbv)  .  pKB)  -.PV)]  ,(r*  .  r2}  w 

b.  Radial  Motion.  To  arrive  at  the  radial  equation  of  motion  it 
is  useful  to  consider  a  volume  with  the  full  ceil  dimensions  Az  and  Ar 
in  the  z  and  r  direction,  but  with  a  small  angular  dimension  A0 
instead  of  2 n,  see  Fig.  2.  Then  one  can  compute  the  radial  motion  from 
Eq.  9  since  that  motion  corresponds  to  a  fixed  direction  in  space.  For 
use  i  n  Eq.  9>  then,  u^  =  u,  and  n^  =  1,  +  1,  A0/2,  0,  0  correspond  to  the 
radial  components  of  the  unit  normals  to  the  inner,  and  outer  surfaces, 
the  two  side  surfaces  and  the  top  and  bottom,  respectively. 


Fig.  2— Element  of  volume  for  discussion  of  radial  motion 
The  pressure-integral  contribution  in  Eq.  9  becomes 

~  A-^mu)  *  Pj  r±  A0  Az  At  -  Pr  r£  A0  Az  At  +  Pg  Ar  A®  A0  At 

where  Pf,  P  and  P  are  the  left,  right  and  side  face  pressures  res- 
x  r  8 

pectively  and  are  defined  in  terms  of  time  N  cell  pressures  as  follows 
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In  terms  of  these  cell  pressures,  and  using  Ar  =  r^  -  r^,  the  equation 
for  radial  motion  becomes 

ijM  -  i  2n  A*  At 

c.  Energy  Equation,  Eg.  10.  In  computing  the  Phase  1  energy  e- 

qiEtion  we  depart  somewhat  from  the  standard  procedure  whiqh  is  used  to 

compute  the  other  integrals  on  the  right  in  Eqs.  7-10.  Specifically, 

2  2 

rather  than  compute  the  gain  in  the  total  cell  energy  ,m(l  +  ^(u  +  v  )) 
by  integrating  the  pressure  forces  on  the  cell  surface,  the  original  e- 
qjation  is  reformulated  to  give  a  volume  integral  expression  for  the  gain 
in  internal  energy,  the  kinetic  energy  increment  having  already  been  fixed 
by  the  preceding  calculation  of  the  velocity  increments.  The  resulting 
equation  is  preferred  because  it  is  a  simpler  expression  for  I  and  be¬ 
cause  it  has  been  used  in  most  of  our  calculations  and  has  led  to  satisfactory 
results.  Starting  with  the  Phase  1  part  of  Eq.  10, 

A-JmE)  =  -  J 8  Pu^ds  ■ 

and  since  cell  mass  m  is  constant  in  Phase  1  md  E  ■  I  +  $  u^u. ,  one 
can  write 


A1(ml)  +  u±  A-^ffluJ  =  -  Js  Pu^ds 


-  Ir  <rui>4T 


OU 

-  r  P  dT  -  J  u,  §§-  &t 

Jt  Ox.^  *>t  i  dxi 

du, 

-  P  f  r— — •  dT  •  U  f  dT 

Ox,  i  *»t  ox7 


where  P  and  in  the  last  line  are  average  values  over  the  cell  volume. 
But  from  the  Phase  1  momentum  equation 


Ai(mui)  *  -  J8  Pr^ds 

A  (»u  )  .  .  J  g-  4t 


so  that  the  last  terms  on  each  side  of  the  energy  equation  cancel 
the  resulting  expression  for  the  gain  in  internal  energy  is 

du. 

^(ml)  -  -P  X,  ^  4t 

xhe  difference  form  is 


SlJ  «(r2  -  rx)  As 

where  (&ui/3x1)N+^  Is  obtained  by  averaging  the  tine  K  (pre- Phase  1)  and 
time  (N+l)  values,  i.e.. 


mhl 


•ht 


P"(K)( 


and,  for  example,  the  time  N  value  of  this  divergence  is  determined  from 
the  calculation 


1  ±  (ru)  ♦  |r 

r  dr  v  7  dz 


.  i  (Tr  U"(KR)  ~  ri  uN(i{l>,\  .  v*(m)  -  v”(kb) 
\hxlJ  rc  \  2dr  J  2lz 

in  which  rQ  *  (rx  +  r2)/2,  ry  =  r2  +  hr/2,  r£.  -  ^  -  hr/2  are  the  radii 
of  cell  centers  for  the  central,  right  ami  left  cells  respectively. 


This  completes  the  Phase  1  calculation,  for  a  typical  cell,  of  the 
momentum  and  internal  energy  increments  due  to  the  pressure  forces  on  the 
cell.  The  two  momentum  increments  are  vised  to  determine  the  velocity 
increments  hy  dividing  by  cell  mass  and  one  has,  at  the  end  of  Phase  1, 
values  of  l(K),  u(K),  v(k)  as  updated  by  Phase  1. 
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3.3.  PHASE  3  -  THE  EFFECT  OF  SHEAR  STRESS 


3.3.1.  Continuity  Equation,  Eg.  8 
No  contribution  in  phase  3. 

3.3.2.  Equation  of  Motion,  Eg.  9 

*•  a.  Axial  Motion.  The  Phase  3  contribution  from  Eq.  9  is 

43(Bu3)  .  at  ;a  e13.yu 

and  for  the  axial  motion 


U3  s  V 

Here  s^ni  is  the  axial  (3)  component  of  the  total  stress  on  a  surface. 
For  a  torus  with  rectangular  section,  s^n^  on  the  various  faces  is 

Bt  •  cell  top 
zz 

-  s*  cell  bottom 
zz 


r  1 

s  _  outer  cell  surface 


-  s^  inner  cell  surface 
rz 


where  zz  and  rz  subscripts  denote  normal  and  shear  stresses  at  the  top, 
bottom,  right  and  left  surfaces  which  are  indicated  by  t,  b,  r  and  l 
respectively.  The  equation  for  axial  motion  is  therefore 

A3(mv)  =  (s^z  -  szz)  TT(rg  -  rj)  At  +  s^z  2nrg  Az  At  -  s*z  2^  Az  At 


b.  Radial  Motion.  For  the  radial  motion  consider  the  element  of 
mass  depicted  in  Fig.  2.  Then  in  the  Phase  3  part  of  Eq.  9, 


we  have 


A3(mu1)  «  At  Ja  s^ds 
u^  ■  u 
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The  radial  components  s^n^  of  the  stresses  on  the  various  faces  of  the 
mass  element  are 

szr  top  of  106188  ei-emeri't 


-  szr  bottom 


s  right 
rr 


-  s*  left 
rr 


-  s 


-  s 


iA0 

’ee 

-^A9 

ee 


front 


back 


The  mass  of  the  element  is  mA6/2n  where  m  is  the  cell  mass.  Tht>n 

&3(mu)  =  (szr  -  szr)(r|  -  rj)  A9  At  +  s^.  rg  A0  Az  At 

\ 

-  rx  A6  At  -  8e0  Az  Ar  A0  At 

oj;  multiplying  by  2n/A0,  the  equation  for  radial  motion  is  . 

A3(mu)  =  21  At  [£(szr  -  szr)(r2  -  ^  *2  Az  -  s^  rx  Az  -  s^AZAr] 

V 

3.3.3.  Energy  Equation,  Eq.  10 

In  the  Phase  3  term  from  Eq.  10, 

A3(mE)  =  At  J8  s^u^ds  , 

s^jU^n^  is  the  work  rate  per  unit  surface  area,  which  for  the  various 

faces  of  the  torus,  is 

/  t  t  t  t.  ,,  , 

(s  v  +  s  u  )  cell  top 


,  b  b  b  b» 

-  (b  v  +  s _ u  )  cell  bottom 

'  zz  zr 
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(Sr  vr  +  Sr  ur)  outer  cell  surface 
rz  rr 

-  i  i  i  i. 

(s  v  +  s  u  )  inner  cell  surface 
rz  rr 

where  t,  b,  r,  L  denote  stresses  and  velocities  at  the  top,  bottom,  right 
and  left  cell  faces  respectively.  The  calculation  of  these  stresses  from 
the  material  constitutive  relations  is  discussed  in  Section  5  and  the 
interface  velocities  are  determined  from  a  simple  average  of  the  time  N  . 
velocities  in  the  two  cells,  e.g. 

yt  ,  vM(K)  +  vN(KA) 

2 

ni  =  IJ^(K)  +  uN(KL) 

2 

The  equation  for  the  Phase  3  change  in  total  cell  energy  is 

.  t  -r»  \  /  t  ■ t  tt  bb  b  b\  /  2  2.  . 

A3(mE)  -  (szz  v  +  szr  u  -  szz  v  -  u  )  n(r2  -  rj  At 

+  (s£z  vr  +  s£r  ur )  2n  r2  Az  At  -  (s*z  vl  +  s^,  u*)  2*  r^  Az  At 

This  completes  the  Phase  3  calculation,  for  a  typical  cell,  of  the 
momentum  and  total  energy  increments  due  to  the  deviator  stresses  acting  on 
the  cell.  The  two  momentun  increments  are  used  to  determine  the  velocity 
increments  by  dividing  by  cell  mass.  This  fixes  the  Phase  3  change  in  cell 
kinetic  energy  £u^u^  so  that  the  internal  energy  (i)  can  be  calculated 
using  the  updated  value  of  total  energy  E  =  I  +  £u^u^  .  At  the  end  of 
Phase  3  one  therefore  has  new  values  of  l(K),  u(k),  v(k)  as  updated  by 
Phase  3* 

3.4.  PHASE  2  -  THE  EFFECT- OF  TRANSPORT 

The  purpose  of  this  section  is  to  describe  how  the  transport  of  mass, 
momenti^n  and  energy  from  cell  to  cell  is  accounted  for  in  the  code.  This 
is  done  by  calculating  the  integrals  in  the  last  terms  of  Eqs.  8  to  10. 

The  method  is  that  given  previously  as  Ref.  2  and  is  a  continuous  analog 
of  the  transport  which  has  been  discussed  in  relation  to  earlier  PIC 
codes. 
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.  3*4.1.  Continuity  Equation,  Eq»  8 
The  transport  mass  is 

A2(m)  -  -  At  Jg  pu^ds 

and  is  determined  for  each  of  the  cell  faces  from 

N  N  .  .. 

om  *  -  p  At 

where  p  is  the  density  of  the  cell  from  which  the  mass  moves  (donor  cell), 

N 

A^  is  the  area  of  the  face  and  is  an  interpolated  value  of  the  velocity 
component  normal  to  A^  representing  approximately  the  velocity'  at  the  inter¬ 
face  at  the  end  of  the  time  step*  For  example,  considering  cell  K  and  the 
cell  KA  above,  one  has 


The  calculated  transport  masses  are  subtracted  from  the  donor  cell  mass  and 
added  to  the  acceptor  cell  mass.  This  updating  is  done,  however,  after  the 
transport  terms  have  been  calculated,  so  that  all  of  the  transport  terms 
are  computed  using  time  N  quantities. 

3*4.2.  Equation  of  Motion,  Eq.  9 

a.  Axial  Motion.  The  term  in  Eq.  9  for  axial  momentum  transport  is 
AgCrni^)  -  -  At  Js  (piu  u±)  nA  ds  . 

At  each  face  of  the  cell  the  transport  specific  momentum,  u^,  is  taken 
to  be  the  axial  velocity  of  the  cell  from  which  the  mass  moves  (donor  cell, 
index  KD),  i.e. 

u j  -  vN(KD)  . 

Since  the  various  faces  of  the  cell  have  different  donor  cells  it  is 
convenient  to  express  the  momentum  transport  for  each  face 

&(mv)  ■  v^(KD)  8m 


14 


where 


*  „N  N  .  . . 

6m  =  -p  u^  A ^  At 

is  the  mass  which  is  transported  across  the  interface,  as  given  in  the  pre¬ 
ceding  mass  transport  calculation.  Note  that  &(mv)  and  6m  are  the  momentum 

© 

and  mass  transports  in  either  the  axial  or  radial  directions,  depending  on 
which  type  of  interface  is  being  computed. 

b.  Radial  Motion.  Again  we  have  from  Eq.  9 

^(rnv^)  *  -At  J  (p^  u^)  n±  ds 

and,  by  analogy  with  the  axial  case,  the  equation  for  the  transport  of 
radial  motion  across  an  interface  is 

6(mu)  =  uN(KD)  6m 

where  uN(KD)  is  the  time  N  velocity  of  the  donor  cell  and  6m  *  -pN  u^A^  At 
is  the  mass  which  is  transported  across  the  interface  in  question,  as  com¬ 
puted  in  3*4.1  above. 

3*4.3  Energy  Equation,  Sq.  (10) 

The  Eq.  (l)  expression  for  the  transport  of  energy  is 

A^mE)  **  -At  J'  (pu^E)  ni  ds  . 


To  evaluate  this  integral,  the  transported  specific  energy  is  taken  to' be 
that  of  the  donor  cell  KD,  i.e. 


ED  s  i^(kd)  +  \ 


L.  _  J 

and  the  total  energy  which  is  transported  across  a  given  interface  is 
therefore  the  product  of  this  specific  energy  and  the  associated  transport 
mass  which  was  computed  above  in  3»4.1., 


A^CmE)  *  Ep  6m 

Once  the  mass,  momentum  and  energy  transports  are  known  for  all.  faces  of 
the  cell,  the  cell  quantities  can  be  updated  for  these  effects.  New  cell 
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velocities  are  determined  by  dividing  the  new  momenta  by  updated  cell  mass. 
This  fixes  nhe  new  kinetic  energy  so  that  the  internal  energy  can  be  cal¬ 
culated  from  the  known  total  energy.  At  the  end  of  Phase  2  one  has  final 
values  of  u(K),  v(k),  I(k),  M(K)  ,  the  cell  velocities,  specific  internal 
energy  and  mass. 

3*5*  CELL  PRESSURES  AND  TIME  STEP  FOR  NEXT  CYCLE 

New  cell  pressures,  to  be  used  in  the  next  time  step,  are  computed 
from  an  equation  of  state  giving  cell  pressure  as  a  function  of  density  and 
specific  internal  energy  within  the  cell.  Such  an  equation  of  state  is  the 
subject  of  Section  4. 

For  stability  of  the  solution  the  new  time- step  is  taken  to  be  less 
than  the  time  for  either  a  sound  wave  or  a  mass  element  to  cross  a  cell. 

In  most  problems  the  actual  time  step  is  some  factor,  0.4  to  0.7,  of  this 
minimum  transit  time. 

3.6.  SPECIAL  FEATURES  OF  THE  COMPUTATION 
3*6.1.  Grid  Boundaries  end  Axis 

\ 

Additional  specification  is  required  for  cells  which  border  the  grid 
boundaries  or  axis  of  symmetry,  because  necessary  quantities  are  not  defined 
for  neighbor  cells  which  would  be  outside  the  grid. 

The  code  provides  for  a  transmittive  boundary  at  the  top  and  right,  an 
option  of  a  transmittive  or  reflective  boundary  at  the  bottom  and  a 
reflective  boundary  at  the  axis.  Boundary  conditions  for  border  cells  are 
then  derived  by  assuming  (fictional)  neighbors  cells  outside  the  grid.  For 
transmittive  boundaries  the  flow  variables  are  the  same  in  the  fictional 
cell  as  in  the  border  cell,  and  for  reflective  boundaries  the  state  is 
assumed  to  be  the  same  except  that  the  velocity  component  normal  to  the 
boundary  has  opposite  sign. 

3*6.2.  Free  Surface  Motion 

Unless  additional  provisions  were  made,  the  preceding  method  would 
result  in  a  diffusion  of  the  free  surface  over  several  cell  dimensions.,  To 
avoid  this,  a  revised  calculation  is  made  to  determine  transport  masses  at 
free  surfaces,  as  described  in  the  next  two  paragraphs. 
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Consider,  for  example,  the  upward  motion  of  a  free  surface,  from  the 
cell  K  which  contains  it  to  the  empty  cell  KA  above.  If  the  density  of  the 
average  cell  K  is  less  than  the  normal  and  if  the  specific  internal  energy 
of  both  cell  K  and  the  cell  below  are  less  than  that  required  to  bring  the 
material  to  its  vaporization  point,  then  the  transport  mass  from  K  to  KA  is 
set  to  zero.  Otherwise,  the  mass  transported  from  K  to  KA  is  given  by 
6ro  -  p(K}  v(K)  A  At  where  A  is  the  area  of  the  interface.  This  procedure 
amounts  to  requiring  that  cell  K  be  filled  before  material  is  transported 
to  KA,  if  the  material  is  condensed.  The  same  procedure  is  used  for  the 
motion  of  a  free  surface  into  empty  cells  to  the  right,  left  or  downward. 

A  similar  calculation  is  made  for  a  receding  free  surface,  i.e.,  one 
which  empties  the  cells  from  which  it  passes.  Consider  the  case  similar  to 
the  one  above,  except  that  the  motion  is  downward,  so  that  one  is  concerned 
with  evaluating  the  transport  mass  from  K  to  KB.  If  either  K  or  KB  has 
specific  internal  energy  greater  than  that  required  to  bring  the  material 
to  its  vaporization  point,  the  normal  mass  transport  calculation  is  made, 
using  the  donor  cell  (K)  density  and  an  interpolated  velocity.  But  if  both 
cells  are  cold,  the  mass  transported  is  5m  *  p(KB)  v(KB)  A  At.  Here  the 
use  of  the  generally  higher  density  p(KB)  causes  a  greater  transport  mass, 
making  it  possible  to  evacuate  free  surface  cells.  In  the  unmodified 
method  some  residual  mass  would  always  be  left  in  an  evacuating  cell.  With 
the  existing  method  the  cell  may  overempty,  but  this  is  immediately 
corrected  should  the  cell  mass  become  negative. 

3.6.3.  Rezone 

A  rezone  subroutine  can  be  triggered  by  mass  motion  out  of  the  right 
or  top  grid  boundaries .  Four  cells  are  combined  into  one,  starting  at  the 
lower  left  corner  (z  =  0,  r  «  0)  of  the  grid  and  new  cells  are  added  at  the 
top  and  right.  The  total  number  of  cells  is  held  constant,  as  well  as 
their  shape. 

3.6.4.  Tracer  Points 

A  system  of  "tracer  points"  can  be  deployed  in  the"  initial  material 
configuration  and  these  points  are  transported  with  local  material  velocity 
during  the  computation.  While  these  points  do  not  enter  the  actual  solution 
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to  the  flow  equations,  their  dispositions  at  various  times  can  give  quick 
and  valuable  insight  into  the  meaning  of  the  solution. 

3.6.5.  Variable  Zoning 

Although  the  discussion  has  proceeded  with  no  special  reference  to 
variable  zoning,  the  current  version  of  OIL-RPM  reported  in  PART  THREE  is 
able  to  compute  problems  in  which  the  zone  size  is  variable.  If  constant 
zones  are  not  specified,  the  cell  sizes  must  be  supplied  as  a  part  of  the 
SETUP  deck.  It  is  possible  to  call  for  variable  cell  size  in  one  direction 
and  constant  in  the  other. 

The  zones  must  be  selected  so  that  the  discontinuities  in  size  between 
adjacent  cells  are  not  excessive  and  the  aspect  ratio  of  the  cells  are  not 
too  extreme.  It  was  found  that  an  initially  spherical  explosion  loses  its 
symmetry  in  the  course  of  a  calculation  when  the  ratio  of  the  cell  edges  is 
greater  than  four  to  one,  and  that  satisfactory  symmetry  is  obtained  when 
the  ratio  is  two  to  one.  Cell  size  variations  are  usually  held  to  less 
than  10  percent  between  adjacent  cells.  This  resulted  in  good  accuracy  in 
a  test  problem. 
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4.  EQUATION  OF  STATE  FOR  THE  CALCULATION  OF  PRESSURE 


Since  the  flow  equations  contain  the  pressure,  p,  the  density,  p,  and 

the  specific  internal  energy,  I,  it  is  natural  to  make  use  of  an  equation 

of  state  relating  these  particular  variables.  A  description  of  such  an 

equation  p(l,p)  developed  for  hypervelocity  impact  calculations  is  given  by 

(9) 

Tillotson.  His  formulation  fits  Thomas-Fermi  data  at  high  pressures 
(10  megabars  and  above),  experimental  shock- wave  measurements  at  more 
moderate  shock  pressures,  and  at  low  density  and  high  energy  it  describes 
a  material  which  behaves  in  the  limit  as  an  ideal  gas.  This  is  done  by 
means  of  an  equation  of  the  form 


At  high  internal  energy  and  low  density  the  equation  of  state  has  the  form 


P  ■  PG(I>P)  »  alp  + 


- -  +  An  e  H'T1 

“2  1 

E0r 


-b(I  -  ’  )  I  -«<4  - 1) 


*  'I 


The  equation  of  state  described  by  Tillotson  is  used  in  a  slightly 
modified  form  in  the  current  program.  The  regime  in  which  material  is 
neither  vaporized  fully  (I  <  I "  )  nor  completely  condensed  (I  >  I#)  is  com¬ 
puted  by  weighting  the  two  expressions  above.  Specifically,  if  I*  <  I  <  I/# 
then 


Here  l'  is  the  energy  of  the  material  which  just  brings  it  to  the  vapor 
temperature  and  l"  includes  the  additional  energy  to  complete  the 
transition  to  the  vapor  state. 
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This  formulation  allows  materia]  to  condense  from  the  vapor  state 
into  the  "dust"  state  in  a  continuous  manner.  By  "dust",  here,  is  meant  a 
state  where  the  internal  energy  is  inadequate  to  bring  material  to  the 
vapor  temperature  but  so  expanded  that  it  cannot  support  pressure. 

In  the  notation  above  I  is  used  to  denote  specific  internal  energy, 

(9) 

as  in  the  other  portions  of  this  report.  TiJJotson, however,  uses  E, 
end  in  portions  of  the  program  documented  in  PART  THREE,  "E"  may  denote 
specific  internal  energy. 

Material  is  allowed  to  support  tensions  (negative  pressures)  in  the 
current  version  of  the  program.  The  maximum  expansion  is  specified  by 
AMDM.  If  the  material  is  cold  (I  <  I/)  and  T)  <  AMDM  the  pressure  is  set  to 
zero.  Then  the  material  is  thought  of  as  cracked,  with  the  cracks  smaller 
in  *  , : e  than  a  cell,  or  as  dust,  depending  on  the  density  of  the  material. 
AMDM  is  typically  a  number  like  0.97*  For  small  expansions  the  expressions 
above  for  pressure  reduce  approximately  to 

p  «  GQ  Ip  +  Am.  , 


where  GQ  is  Gruneisen's  ratio.  \ 

Thus,  for  aluminum,  assuming  the  first  term  is  negligible,  we  would 
find  p  »  -22.5  kiiobars  using  the  value  A  *  O.75  megabars  given  by 
Ti  Hots  on. 

A  table  of  constants  appropriate  for  representing  a  variety  of  metals 


is  given  in  PART  THREE.  The  information  for  metals  is  taken  from  the 
Tillotson  report  and  for  geologic  materials  it  is  taken  from  a  report  by 
Allen. ^ 10 ^ 


Some  general  remarks  can  be  made  regarding  the  above  equation  of  state, 
a.  The  pressure  is  continuous  in  p  and  I  in  the  various  regions  and 
across  transitions  between  the  regions.* 


*  We  are  indebted  to  Caroti  of  General  Electric  for  pointing  out  to  us  an 
error  in  an  earlier  version.  Specifically,  isentropes  computed  from  the 
equation  of  state  given  in  Ref.  9  can  be  discontinuous  in  pressure  at  the 
condensed-vapor  transition.  The  correction  of  this  error  has  lead  to 
significantly  smoother  solutions  to  several  applications  of  RPM. 
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b.  No  equation  of  state  formulation  is  ideally  suited  to  all  appli¬ 
cations.  The  present  one  is  especially  suited  to  situations  in  which 
the  material  experiences  a  strong  shock  (say  tens  of  megahars)  and 
subsequently  expands.  Examples  are  the  hypervelocity  impact  problem 
and  certain  ground  shock  applications.  For  quite  different  appli¬ 
cations,  other  equation  of  state  formulations  should  be  considered. 

For  example,  problems  involving  cold  highly  compressed  states 
(p/Pq  >  3,  not  achieved  in  a  single  shock)  should  be  computed  using 
an  equation  of  state  which  is  formulated  with  such  states  in  mind. 

Also  simpler  equations  of  state  may  be  desired  in  some  cases.  An 
example  is  ordinary  shock  wave  hydrodynamics  at  pressures  less  than  a 
megabar,  where  several  simpler  forms  have  been  widely  used  in  past 
years. 

It  is  expected  that  other  losers  of  OIL-RPM  may  want  to  Include  addi¬ 
tional  equations  of  states.  '  No  special  difficulties  are  envisioned  in  such 
extensions. 
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5.  CONSTITUTIVE  EQUATIONS  FOR  THE  CALCULATION  OF  DEVIATOR  STRESS 

In  Section  3.3  the  effect  of  shear  stresses  on  the  calculation  of 
velocity  and  internal  energy  increments  (through  the  conservation  equations) 
was  discussed.  The  total  stress  was  broken  up  into  an  isotropic  component 
(pressure)  and  a  deviator  stress.  The  calculation  of  pressures  from  density 
and  internal  energy  was  then  discussed  in  Section  4.  The  finite  difference 
technique  for  calculating  the  deviator  stresses  from  the  velocity  field 
will  be  the  subject  of  this  section.  It  is  based  on  the  rigid-plastic 
model  relating  deviator  stress  and  deviator  strain  rate  discussed  in  PART 
ONE.  The  motivation  for  selecting  the  rigid- plastic  model  (primarily  its 
simplicity)  among  the  possible  bases  for  a  constitutive  equation  was  also 
described  there. 

5.1.  FINITE  DIFFERENCE  APPROXIMATIONS 

Deviator  stresses  act  on  each  face  of  a  cell  to  accelerate  the  mass  it 
contains,  as  sketched  in  Fig.  3*  A  tabulation  of  these  stresses  is  gi  j. 
below,  and  in  the  table  the  relation  of  the  current  notation  to  that  used 
in  Section  3,  is  also  indicated.  This  current  notation  is  consistent  with 
the  FORTRAN  symbols  used  in  the  program  in  PART  THREE  of  this  rep  "t. 


TABLE  5.1 

NOTATION  FOR  STRESSES  ACTING  ON  A  CELL 


Area 

Normal 

Tangent 

Hoop 

Top 

SNT  « 

St 

zz 

STT  * 

Szr  - 

— 

Right 

SNR  » 

sL 

rr 

STR  « 

4 

— 

Left 

SNL  - 

s* 

rr 

STL  = 

4 

- — 

Bottom 

SNB  ** 

Sb 

zz 

STB  - 

4 

— 

Cross 

Section 

— 

• 

J 

1 _ 

Hoop 

Fig.  3 — Deviator  stresses  used  in  the  Phase  3  calculation 


The  notation  in  the  table  is  related  to  the  more  formal  subscript 
notation  of  PART  ONE  by  the  equations 

S11  =  Srr  '  S13  =  Srz  '  S22  =  SG6  '  S33  *  Szz  ' 

The  required  deviator  stresses  are  computed  from  the  strain  rates  by 
means  of  the  relation 

ht-'W 
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indicated  in  PART  ONE,  where  Y  is  the  yield  strength  in  simple  shear. 


eij  “  eij  "  95ij  ' 


9  *  cii^3 


and 


M  ■  eu  eu 


<•«> 


ur 

0 


-a W 


The  strain- rate  tensor  in  cylindrical  coordinates  is  given  by  the  matrix 

0 

u/r  0 

0  v 


?<VTr> 


where  the  subscripts  r  and  z  denote  differentiation  with  respect  to  the 
subscripted  variable.  Expressions  for  6  and  W  are  obtained  by  substitution 
of  the  matrix  elements  into  the  expressions  above. 

0  -  -~(u  +  -  +  v  ) 

3  r  r  z 

W-u  2  +  v  2  +  (u/r)2  +  -i(u  +  v  )2  -  302  • 
r  z  '  2'  z  r 


The  velocity  gradients  vised  in  calculating  the  stresses  are  deter¬ 
mined  in  such  a  way  that  the  velocities  of  the  cells  adjoining  each  inter¬ 
face  enter  into  the  calculation  in  a  symmetrical  manner.  Details  of  the 

differencing  are  indicated  in  Table  5.2.  The  cell  size  is  DRV(l)  in  the 
radial  direction,  DZV(j)  in  the  axial  direction, 
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TABLE  5.2 

CALCULATION  OF  VELOCITY  GRADIENTS 


Top- centered  Variables 


du  _  u(KA)  -  u(K) 
dz  "  DZV(j) 


dv  V(KA)  -  v(K) 
dz  =  DZV(j) 


au  |(u(KR)  +  u(KAR))  -  -|(u(KL)  +  u(KAL)) 
dr" =  2DRV(I) 

|(v(KR)  +  v(KAR))  -  |(v(KL)  +  v(KAL) ) 
dr  =  2DRV ( I ) 


u  u(KA)  +  u(k) 
r  *  2R( I ) 

DZV(J)  =  -|(DZ(j+l)  +  DZ(j)),  2DRV ( I )  -  -§DR(l-l)  +  DR(l)  +  -|dr(I+1) 


Right- centered  Variables 


du 

■|(u(KA)  +  u(KAR) )  -  -|(u(KBL)  +  u(KBR)) 

du 

dz 

EZV(J) 

dr 

dv~ 

-|(v(KA)  +  v(KAR))  -  -|(v(KB)  +  v(KBR)) 

dv 

dz 

DZV(j)  “ 

dr 

u  u(K)  +  u(KR) 
r  =  R(l )  +  R(l+1) 

DRV(I)  =  -|(DR(I)  +  DR(l+l) ), 

DZV(J)  » 

„  y.(KR)  -  u(K) 
DRV(l) 

.  y(KR),  v(Kl 
DRV(l) 


DZ(J) 


Cell- centered  Variables 


du  _  u(KA)  -  u(KB) 
dz  ”  2 DZV(j) 

dv  a  v(KA)  -  v(KB) 
dz  "  2DZV(J) 


du  u(KR)  -  u(KL) 
dr  "  2DRV ( I ) 

dv  =  v(KR)  -  v(KL) 
dr  =  2DRV ( I ) 


u_  u(K) 
r  =  R(I) 

2  DZV  =  -|dZ(J+1)  +  DZ(J)  +  -|dZ(J-1) 
2DRV  -  -|dr(I+1)  +  DR(I)  +  ^DR(l-l) 


and,  as  suggested  by  the  notation,  these  may  vary  from  cell  to  cell. 
Expressions  for  DRV(l)  and  DZV(j)  in  terms  of  cell  dimensions  DR(l)  and 
DZ(j)  are  included  in  the  table.  DR(l)  and  DZ(  J)  are  abbreviated  as  DR 
and  DZ  in  some  of  the  following  equations  for  brevity  of  notation. 

Values  for  the  left  and  bottom  deviator  stresses  do  not  have  to  be 
determined  separately  since  they  are  calculated  as  the  right  and  top 
deviator  stresses  for  the  cells  to  the  left  and  below  respectively.  In 
the  finite  difference  notation  described  in  the  preceeding  paragraphs  the 
equations  of  Section  3*3  become 

A^mv)  »  2nAt  [ (SNT  -  SNB)  •  R*  DR  +  (STR.R**  -  STL*R*)  •  DZ] 

and 

A^nu)  »  2nAt  [Rr  •  DR(STT-ST'1)  +  DZ(SNR*Rr  -  SNL*R*  )  -  HOOP.DR'DZ} 

vi  § 

where  R  is  the  coordinate  of  the  cell  center  and  R  and  R  are  coordinates 
of  the  right  and  left  edges: 

Rr  «  R  +  -|dR,  R*  -  R  -  -|dr  . 

To  calculate  the  change  in  internal  energy  in  a  cell  it  is  noted  that 
the  change  in  its  total  energy  is  determined  by  the  stresses  acting  on  its 
boundary,  and  the  change  in  kinetic  energy  is  determined  by  the  change  in 
velocity,  which  in  turn  is  obtained  from  the  momentum  equations.  The 
difference  of  the  changes  in  total  and  kinetic  energy  in  the  change  in  inter¬ 
nal  energy.  The  change  in  internal  energy  cf  a  cell  due  to  pressure  is 
accounted  for  in  Phase  1,  and  the  change  due  to  the  deviator  stresses,  which 
may  be  thought  of  as  the  plastic  work,,  is  accounted  for  in  Phase  3*  The 
Phase  3  change  in  internal  energy  is  computed  in  several  parts,  viz, 

h3(ml)  =  (A+B)  At  -  [(u+-|a3u)  A3u  +  (v+^v)  AjV] 

where 

A  -  2ttR-DR  ||(u(k)  +  u(KA))  •  STT  +  |(v(k)  +  v(KA)).  SNT 
-  |(u(K)  +  u(KB))*STB  -  !(v(K)  +  v(KB))  •  SNb| 

2£ 


ijM. 


and 

B  -  2n  DZ  |-|Rr  .  (u(KR)  +  u(K))*SNR  +  *|  Rr  (v(KR)  +  v(K))‘STR 

-  -|r£(u(K)  +  u(KL) )•  SNL  -  ^R£(v(k)  +  v(KL))*  STl| 

The  plastic  work  done  in  a  single  cycle  in  a  single  cell  is  It  is 

positive  in  principle,  though  not  necessarily  in  the  calculation,  as  a 
result  of  the  finite  time  step.  The  sum  of  the  plastic  work  done  on  each 
cell  in  each  cycle  is  accumulated  into  a  grand  total  which  is  printed  out 
at  each  of  the  edit  times.  This  is  the  "plastic  work"  exhibited  in  the 
edit  prints. 

It  is  possible  for  the  change  in  velocity  in  a  single  cycle  due  to 
plastic  stresses  to  be  so  large  that  the  velocity  profile  changes  curvature 
This  difficulty  can  be  largely  corrected  by  subcycling.  Specifically,  the 
time  step  for  Phase  3>  is  divided  into  N  steps,  where 

A1  “  N  '  V  "  (N_1)  N(N+i)'  '  V  *  N(N+l) 

The  number  of  time  steps  (CYCPH3)  is  left  to  the  discretion  of  the  user, 
but  experience  has  shown  that  4  is  a  good  choice. 

The  problems  that  the  OIL-RPM  code  has  been  called  upon  to  calculate 
have  involved  shocks  diverging  out  from  either  an  impact  or  an  explosion. 

In  either  case  a  very  significant  improvement  is  obtained  by  preventing 
the  Phase  3  calculation  from  affecting  the  flow  outside  the  shock.  To  do 
this  the  shock  has  to  be  located,  at  least  in  an  approximate  manner.  This 
is  done  by  seeking  out  the  location  of  the  first  pressure  maximum  in  each 
column  of  cells.  The  result  is  used  to  update  the  JPM(l)  array,  which 
gives  the  J  index  of  the  cell  with  the  Pressure  Maximum  in  column  I.  The 
updating  of  JFM(l)  is  done  in  the  subroutine  CDT  which  determines  the  time 
step.  The  JPM(l)  have  to  be  halved  when  rezoning,  and  that  is  done  in  the 
subroutine  "REZONE." 
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5.2.  VARIATIONS  II :  YIELD  STRENGTH 

Experiments  (ll)(l2^ave  shown  that  yield  strength  depends  on  pressure 
and  temperature .  Although  an  appropriate  set  of  constitutive  equations 
fully  describing  these  effects  is  not  available,  the  calculation  of  Phase  3 
attempts  to  account  for  them  by  allowing  the  flow  stress,  Y,  to  vary  ac¬ 
cording  to  the  relation 

Y  -  (Yq  +  Yxn  +  Y^2){1-I/e0) 

where 


I  is  the  specific  internal  energy  and  EQ  is  the  energy  which  raises  the 
material  to  the  melting  point.  If  either  factor  is  negative  the  yield 
strength  is  set  to  zero.  This  formula  is  not  intended  as  an  accurate 
representation  of  the  physical  behavior  of  materials,  but  was  incl'.ided  in 
the  program  so  that  the  importance  of  yield  strength  variations  could  be 
estimated. 

Although  no  formal  report  has  been  written,  it  seems  appropriate  to 

* 

mention  here  a  calculation  in  which  the  second  factor  in  the  expression 
above  for  Y  was  set  equal  to  unity.  The  calculation  otherwise  duplicated 
the  standard  impact  crater  problem  reported  in  Ref.  14.  The  effect  was  to 
decrease  the  crater  depth  at  lU  y,sec  from  1.16  to  1.11  cm.  At  that  time 
the  rate  of  change  of  depth  is  down  to  10^  cm/sec  and  thus  its  growth  is 
essentially  terminated.  In  a  second  comparison  it  was  found  that  setting 
the  first  factor  to  unity  made  the  crater  depth  at  14  psec  1.18  cm  rather 
than  the  1.16  cm  computed  in  the  standard  problem. 

'"To  summarize,  the  effect  of  yield  strength  variations  can  be  accounted 
for  in  an  approximate  manner  by  choosing  appropriate  values  for  Y^,  Yg  and 
Eq,  but  this  does  not  seem  to  influence  crater  size  in  the  first  approcti- 
mation.  EQ  can  be  found  from  handbooks  as  the  internal  energy  at  the 
point  where  the  yield  strength  goes  to  zero,  and  this  can  be  taken  roughly 
as  the  internal  energy  at  melting  in  the  absence  of  detailed  informat-ion. 

Y1  and  Yg  can  be  estimated  from  the  relation,  Y  -  YQ  +  a?,  which  is  valid 
at  low  pressures.  The  value  of  Of  is  about  0.07  for  metals  and  near  unity 
for  geological  materials. 


28 


$  o 


REFERENCES 


1 


1.  Rich,  M.,  "A  Method  for  Eulerian  Fluid  Dynamics,"  Los  Alamos 
Scientific  Laboratory  Report  LAMS-2826  (1963) . 

2.  Johnson,  W.E.,  "OIL:  A  Continuous  Two-Dimensional  Eulerian 
Hydrodynamic  Code,"  General  Atomic  Report  GAMD-5580  (1965). 

3.  Gentry,  R.A.,  R.E.  Martin,  B.J.  Daly,  "An  Eulerian  Differencing 
Method  for  Unsteady  Compressible  Flow  Problems,"  J.  Computational 
Physics  1,  (1966) . 

4.  Wilkins,  Mark  L.,  "Calculation  of  Elastic-Plastic  Flow,"  Methods 
in  Computational  Physics,  Vol.  3>  211,  Academic  Press  (1964). 

5.  Maenchen,  G.,  and  S.  Sack,  "The  Tensor  Code,"  Methods  in  Com¬ 
putational  Physics,  Vol.  3>  l8l,  Academic  Press  (1964). 

6.  Dienes,  J.K. ,  W.E.  Johnson,  J.M.  Walsh,  "Annual  Status  Report  on 
the  Theory  of  Hypervelocity  Impact,"  General  Atomic  Report 

GA-6509,  (1965). 

7.  Evans,  M.W.,  F.H.  Harlow,  "The  Particle  in  Cell  Method  for  Hydro- 
dynamic  Calculations , "  Los  Alamos  Scientific  Laboratory  Report 
LA-2139  (1957). 

8.  Harlow,  F.H.,  "Two  Dimensional  Hydrodynamic  Calculations,"  Los 
Alamos  Scientific  Laboratory  Report  LA-2301  (1959)* 

9 .  Tillotson,  J  - "Metallic  Equations  of  State  for  Hypervelocity 
Impact,"  Gen  ral  Atomic  Report  GA-3216,  July  1962. 

10.  Allen,  R.T.,  "Equation  of  State  of  Rocks  and  Minerals,"  General 
Atomic  Report  GAMD-7834,  March  1967. 

11.  Bridgeman,  P.W.,  "Studies  in  Large  Plastic  Flow  and  Fracture," 
Metallurgy  and  Metallurgical  Engineering  Series,  McGraw-Hill,  1952. 

12.  Zhurkov,  S.N.  and  T.P.  Sanfirova,  "A  Study  of  the  Time  and  Temperature 
Dependences  of  Mechanical  Strength,"  Soviet  Physics  Solid  State, 

Vol.  2,  No.  6,  December  i960. 

13.  Dienes,  J.K.  and  M.W.  Evans,  "Cratering  Calculations  with  a  Hydro- 
dynamic  Strength  Code,"  General  Atomic  Report  GAMD-7369, 

September  1966. 


29 


